Out-of-equilibrium processes in suspensions of oppositely charged colloids: 
liquid-to-crystal nucleation and gel formation 
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We study the kinetics of the liquid-to-crystal transformation and of gel formation in colloidal 
suspensions of oppositely charged particles. We analyse, by means of both computer simulations 
and experiments, the evolution of a fluid quenched to a state point of the phase diagram where 
the most stable state is either a homogeneous crystalline solid or a solid phase in contact with a 
dilute gas. On the one hand, at high temperatures and high packing fractions, close to an ordered- 
solid/disordered-solid coexistence line, we find that the fluid-to-crystal pathway does not follow the 
minimum free energy route. On the other hand, a quench to a state point far from the ordered- 
crystal/disordered-crystal coexistence border is followed by a fluid-to-solid transition through the 
minimum free energy pathway. At low temperatures and packing fractions we observe that the 
system undergoes a gas-liquid spinodal decomposition that, at some point, arrests giving rise to 
a gel-like structure. Both our simulations and experiments suggest that increasing the interaction 
range favors crystallization over vitrification in gel-like structures. 

PACS numbers: 



Colloidal suspensions are excellent model systems to 
study condensed matter physics. The reason is that col- 
loids can be observed directly with current experimen- 
tal techniques. In addition, colloids can form structures, 
such as colloidal crystals d] or gels 0, that have potential 
applications, for instance, in photonics or in cosmetics. 
It is therefore essential to gain a better understanding 
in the physics governing the spontaneous transformation 
of a disordered fluid phase into a self-assembled colloidal 
structure. Here, we study how oppositely charged col- 
loids assemble to form either crystals or gels. In partic- 
ular, we are interested in the mechanism by which the 
fluid transforms into either phase. In our work, we make 
use of computer simulations to study the fluid-to-crystal 
transition, whereas the gel formation has been studied 
both numerically and experimentally. 

Quite recently, it has been demonstrated experimen- 
tally that a suspension of oppositely charged particles can 
transform into different crystal phases [lj, and that the 
equilibrium phase diagram of a simple model potential 
for oppositely charged particles exhibits the same crystal 
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structures found in experiments Q. More recently, the 
gas-liquid coexistence line was calculated for the afore- 
mentioned model potential 0] ■ Taking advantage of the 
accurate knowledge of the equilibrium phase diagram we 
study the kinetics of both fluid-solid and fluid-gel transi- 
tions in oppositely charged colloidal suspensions. 

The fluid-to-solid transition is an intriguing problem 
widely studied in colloid physics . The reason is that, 
by understanding the initial steps of crystal growth, one 
can aim to control the size and the structure of the crys- 
tallites eventually formed. The study of crystal nucle- 
ation in suspensions of hard colloidal particles has already 
received quite some attention both experimentally 0, 0] 
and by computer simulations [B], Q . Other colloidal sys- 
tems for which the liquid-to-solid transition mechanism 
has been studied are, for instance, equally char ged par- 
ticles 0, [l(| or binary mixtures of hard spheres fl2| . 
Here, we present the case of oppositely charged parti- 
cles. A special characteristic of the phase diagram of 
oppositely charged spheres is the coexistence between 
ordered and substitutional^ disordered structures on a 
regular 3D lattice [1, U|. This offers a good opportunity 
to study the role of an order-disorder phase transition 
in the liquid-to-solid nucleation pathway. We find that, 
in certain thermodynamic conditions, the formation of a 
disordered phase is kinetically favored, whereas the nu- 
cleation of the ordered phase entails a lower free energy. 

The formation of colloidal gels has also received great 
attention; not only because of the wide use of col- 
loidal gels in cosmetics, drugs or food additives, but also 
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because understanding gelation remains a fundamental 
challenge. The most popular system regardin g co lloidal 



gelation is that of colloid-polymer mixtures [14| . The 
presence of non-adsorbing polymers induces an effective 
attraction between the colloids, which is responsible for 
the gelation to happen 15, [l6j]. For the case of attrac- 
tive particles, the mechanism by which a colloidal fluid 
transforms into a gel has been described as the arrest of 
a gas- liquid spinodal decomposition 0, [13, EH • When a 
low density fluid separates into a very low density gas 
and a high density liquid by spinodal decomposition, a 
percolating pattern of dense liquid forms after the first 
steps of the spinodal phase separation. If the dynamics 
of the particles contained in the network- forming liquid is 
sufficiently slow, the pattern gets arrested, giving rise to 
a gel. Recently, it has been pointed out that mixtures of 
oppositely charged colloids can also form 

gels 0,0, [2l|. 

Hereby, we present experimental evidence confirming this 
statement. Moreover, we make use of computer simula- 
tions to gather evidence that the mechanism of gel for- 
mation for this system is an arrested spinodal decom- 
position, as it is the case for purely attractive particles. 
Finally, we study, both numerically and experimentally, 
the role of the interaction range on the interplay between 
vitrification and crystallization. Our results indicate that 
increasing the interaction range locally favors crystalliza- 
tion over vitrification. 



METHODS AND TECHNICAL DETAILS 
Computer Simulations 

In the present work, we study a symmetric binary 
mixture of equally-sized oppositely-charged particles. 
The number of particles ranges from 686 (for the gel- 
formation simulations) to 8000 (for the Umbrella Sam- 
pling simulations). 1000 particles were used for the spon- 
taneous nucleation and Forward-Flux-Sampling (FFS) 
runs and for some of the gel-formation states. 

In our computer simulations, the screened Coulomb in- 
teraction between two colloids of diameter a and charge 
Ze is approximated by a Yukawa potential plus a repul- 
sive core: 



u{r)/k B T = { ±u* 




; -K(r-o-) 

r/cr 



r < a 

a < r < r c 

r > Tr 



(1) 



where r is the distance between the centre-of-mass of 
the colloids, u* is the energy at contact in k B T units 
and r c is the cut-off radius (set at the distance where 
u(r)/(u*k B T) = 10~ 3 ). The energy in k B T when r = a, 
u*, is equal to Z 2 X B /{{l + ^-) 2 a) 1 where X B = e 2 /e s k B T 
is the Bjerrum length (e s is the dielectric constant of 
the solvent) and n — \/8ir\ B p sa it is the inverse Debye 



screening length for a 1:1 electrolyte (p sa it is the number 
density of added salt ions). To study the liquid-to-solid 
transition, we performed NPT Monte Carlo (MC) sim- 
ulations using Eq. [Tj as the Hamiltonian. In order to 
study the formation of gels, we performed Brownian Dy- 
namics (BD) simulations in the NVT ensemble, where 
we substituted the hard-core interaction of Eq. [I] by a 
steep u* /r 36 repulsion. Our units of energy and length 
are u* and a respectively. Consequently, the time units 
are t* = a 2 / (Dqu*). The time step we used for the inte- 
gration of the Langevin position equation 22] is 7-10~ 6 t*. 

Figure Q] shows a sketch of the equilibrium phase di- 
agram for the model potential described by Eq. [IJ in 
the u*-colloidal-packing-fraction (<j>) plane, calculated for 
recr = 6 d SI- At low and u* , the fluid is the most 
stable phase. At high cf>, there exist several crystal struc- 
tures. Two of them coexist with the fluid: a disordered 
face-centred-cubic solid (disordered fee) at low u* , and 
a CsCl-like structure at high u*. In the disordered- fee 
solid, particles arrange on an fee lattice, but the sign of 
the charge each particle bears does not follow any order 
throughout the crystal. At high u* the fluid-CsCl coex- 
istence opens up and the solid coexists with a very low 
density fluid (gas). Buried underneath the gas-solid co- 
existence, there is a metastable gas-fluid coexistence line, 
which is denoted in the phase diagram by a dashed line. 
In our simulation, we quenched a stable fluid to the state 
points indicated by an asterisk and a letter in Fig.[TJ and 
studied the structural evolution of the fluid in its way 
towards either a solid or a gel. 

In order to observe a rare event such as crystal nu- 
cleation in a simulation of a modest number of parti- 
cles it is necessary to use ad hoc simulation techniques. 
We have used two of them: Umbrella Sampling [23| and 
Forward Flux Sampling 24j. The Umbrella Sampling 
method (US) involves biasing the sampling of an equilib- 
rium MC simulation in order to explore high free energy 
configurations (containing large pre-critical clusters). We 
refer the reader to the review published in [10] for the 
method details. The Forward Flux Sampling technique 
(FFS) uses an order parameter relevant for the transi- 
tion to split the free energy landscape in regions divided 
by hyper-surfaces (or interfaces). In the FFS scheme, 
the calculation of the flux of formation of post-critical 
solid clusters in the liquid (number of post-critical clus- 
ters formed per unit of time and volume) is substituted 
by the flux of formation of small pre-critical solid clusters 
times the probability to form a post-critical cluster once 
one of these small pre-critical clusters has been formed. 
In contrast to the US method, FFS can also be used 
in strongly out-of-equilibrium situations, such as crystal- 
lization under shear (2||. FFS was applied for the first 
time to crystal nucleation to calculate the homogeneous 
crystallization rate in molten NaCl [26| . 

In order to study the liquid-to-solid transition path it 
is necessary to identify the particles belonging to a solid 
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FIG. 1: Phase diagram of equally sized oppositely charged 
particles (taken from Refs. |3j] and 3). The Hamiltonian 
is given by Eq[T] (no = 6). The solid lines indicate the co- 
existence lines, whereas the dashed line corresponds to the 
(metastable) gas-liquid binodal (solid-solid coexistence re- 
gions are too narrow to be visible on the scale of the fig- 
ure) . We have studied the structural evolution of a fluid when 
quenched to the state points indicated by an asterisk accom- 
panied by a letter. The state points corresponding to each 
letter are the following: A:( u = 2, <j> = 0.493), B:(u* = 1, 
4> = 0.553), C:(u* = 1, tj> = 0.526), D:(it* = 6.5, 4> = O.f), 
E:(u* = 8, </> = 0.1), F:(u* = 10, (j> = 0.1), G:(u* 
4> = 0.1), R:(u* = 6.5, <(> = 0.3), 
J:(u* = 10, 4> = 0.3), K:(u* = 15, 
cj> = 0.5), M:(u* = 8, 4> = 0.5), 
0:(u* = 15, = 0.5). 
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cluster. To do that, we use a local bond order parameter 
analysis [27|. In our study we compute the qe complex 



vector for each particle, which is a function of the relative 
orientation of a particle with respect to its neighbours 
(two particles are neighbours if they are closer than 1.33 
a, the first g(r) minimum). The component m of the 



vector associated with the i th particle is given by 



1 ^l,m{j^ij 

% (q?q*?) 1/2 



<ii. 



m 



[-6,6] (2) 



where, for our case, I = 6, iVj is the number of neighbors 
of particle i, Yj m (rjj) is a spherical harmonic function 
whose form depends on / and m and whose value depends 
on the relative orientation of particles i and j (r^ ) . In a 
perfect bcc or fee crystal, all the particles have the same 
environment, and, therefore, the scalar product between 
the vectors of any pair of particles is 1. Two neighbour 
particles are considered to be "connected" if the scalar 
product of their normalised vectors exceeds a threshold 
of 0.7. If a particle has at least 9 connections it is labelled 
as "solid". Any two solid-like particles closer than 1.3 a 
belong to the same cluster. In this way we identify all 
solid clusters present in the metastable fluid phase (the 
size of the biggest one is the order parameter we have 



used to study crystal nucleation) . It is important to note 
that our criterion to identify solid-like particles is blind 
to the type of solid lattice (bcc or fee). This means that 
by biasing the system to grow solid clusters with our 
solid criterion we are not forcing the appearance of any 
particular solid structure. 

To measure the degree of substitutional charge disorder 
within the growing solid clusters we use the following 
charge order parameter: 



i=i 3=1 



(3) 



where N is the number of particles in the solid clus- 
ter, Ni is the number of neighbours of particle i (for this 
purpose, a neighbour is any particle at a centre-of-mass 
distance closer than 1.12c), and a is either +1 or -1 de- 
pending on whether the particle is positively or nega- 
tively charged. Thus, £ is equal to 1 if all particles are 
surrounded by oppositely charged neighbours and -1 if 
all particles are surrounded by equally charged neigh- 
bors. Hence, for a perfectly ordered phase one should 
expect £ = 1, whereas for a random phase £ should be 
0. Of course, when thermal fluctuations are present, £ 
is lower than 1 even if the solid phase is ordered. How- 
ever, on average, £ will always be higher for a substitu- 
tionally ordered phase than for a disordered one. (Note 
that the threshold used to identify neighbouring particles 
here (1.12 a) is different from that used for the orienta- 
tional order parameter described above (1.33 a). Both 
thresholds have been set to maximize the difference in 
order parametre value between either "liquid" and "solid" 
particles or substitutionally ordered and substitutional^ 
disordered clusters). 

When studying crystallization at low and constant den- 
sities (NVT ensemble) we have used a looser criterion to 
detect solid-like particles. The reason is that in those 
conditions there are many particles at the interface be- 
tween a dense phase and a gas phase. Obviously, an in- 
terface particle can not have as many solid "connections" 
as a bulk particle does. For us, a particle is solid-like if 
it is "connected" to at least 50% of its neighbours. Two 
particles are "connected" if the scalar product of their q 6 
vectors is larger than 0.5. The neighbouring distances 
are defined in the same way as they were for the liquid- 
to-solid transition study at high densities (see above). 
We use the NVT ensemble for the study of the system's 
evolution at low densities because we are interested in 
the formation of gel-like structures. These structures are 
formed in (gas-solid) coexistence regions of the phase di- 
agram. Since in &p — T projection of the phase diagram 
there are not coexistence regions -but coexistence lines- 
the iVpTensemble is not suitable if one desires to keep 
the system at coexistence. 
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Experiments 

We used two different batches of polymethylmethacry- 
late (PMMA) particles, which were covalently labelled 
with the fluorescent dyes rhodamine isothiocyanate 
(RITC) and 7-nitrobenzo-2-oxa-l,3-diazol (NBD), and 
sterically stabilized with poly-12-hydroxystearic acid. 
We synthesized these particles by means of dispersion 
polymerization 28]. Both particle types were 2.5^m in 
diameter, with a size polydispersity of 3 %, as determined 
from static light scattering measurements. 

We suspended the particles in a mixture of as received 
cyclohexylbromide (CHB, Fluka) and 26.5% cis-dccalin 
(Sigma- Aldrich) by weight. This mixture nearly matched 
the refractive index and density of the particles and had 
a dielectric constant of 5.G. We made suspensions with 
a 1:1 number ratio of the differently labeled particles 
at an overall volume fraction of (f> 0.20, and either 
with (0.5/zM) or without the charge determining tetra- 
butylammonium bromide salt (TBAB, Sigma- Aldrich) 
[E i<|. In CHB, PMMA particles acquire a positive 
charge. When TBAB salt is added a fraction of the 
ions adsorbs onto the particle, with more of the small 
Br~ ions adsorbing than of the large organic TBA + ions. 
The initial positive charge and the charge sign rever- 
sal point differ for each particle batch. We then made 
gradient samples by filling a glass capillary (Vitrocom) 
partially with the salt-free suspension and partially with 
the 0.5//M TBAB suspension, and allowing it to form 
a macroscopic salt gradient over a couple days' time. 
Suspending the particles in a salt gradient is a conve- 
nient way to quickly explore different experimental con- 
ditions. We studied these samples by means of confocal 
microscopy. 

RESULTS 

Homogeneous crystal nucleation at high 
temperature and constant pressure 

The most widespread view of crystal nucleation is that 
provided by Classical Nucleation Theory (CNT) [lo| ■ Ac- 
cording to CNT, the formation of a crystal nucleus in a 
metastable liquid costs free energy due to the formation 
of an interfacial area A. On the other hand, the fact 
that the solid's (sol) chemical potential (/i) is lower than 
that of the metastable fluid (flu), is the driving force for 
crystal nucleation, and lowers the free-energy of crystal 
formation. The balance between these two terms gives 
rise to a free energy barrier (AG): 

AG = A-i - n\fi sol - fi mf \ (4) 

where 7 is the average liquid-solid surface free en- 
ergy and n is the number of particles contained in the 



largest solid cluster. Because of the nucleation barrier, a 
(metastable) fluid can exist at state points of the phase 
diagram where the solid has a lower chemical poten- 
tial. Given that a metastable fluid is a quasi-equilibrium 
state, equilibrium thermodynamic statistics can be ap- 
plied to study its properties. This allowed Duijneveldt 
and Frenkel to develop the US method to calculate nu- 
cleation free energy barriers [3]]]. The US method re- 
lies on the assumption that pre-critical clusters are in 
quasi-equilibrium with the surrounding liquid. Hence, it 
is possible to use any kind of MC sampling scheme to 
analyse the cluster's properties (structure, shape, criti- 
cal cluster size...) via an US scheme. Here, we use two 
types of MC simulations to study the properties of the 
pre-critical clusters: (i) MC NPT simulations, in which 
volume moves and single particle displacements are used 
to sample the configurational space, and (ii) MC NPT 
simulations with charge swap moves, in which, besides 
the volume and particle moves, we occasionally try (20% 
of the times) to swap the positions of two oppositely 
charged particles. As mentioned before, as long as the 
quasi-equilibrium assumption holds, the properties of the 
clusters should not be affected by the way the system is 
sampled. 

The US sampling scheme is tailored for conditions at 
which there is a high free energy barrier (few tens of 
ksTs) separating the metastable fluid and the solid. If 
the free energy barrier is not too high -nor too small 
(otherwise there would not exist a metastable liquid)- 
one can observe the liquid-to-solid transition without any 
biased simulation technique. This is precisely the case of 
state point A, which is analysed in FigO Fig[2] (a) shows 
the evolution of the internal energy per particle for two 
MC simulations (with and without charge swap moves) 
starting from a state-point-A metastable fluid configu- 
ration. After an initial stage in which the energy fluc- 
tuates around an average value (metastable liquid), the 
liquid transforms into a solid -as the energy drop proves. 
The solid formed in both cases (with and without charge 
swaps) corresponds to the stable phase at the state point 
A: CsCl. Let us now analyse the nucleation pathway 
in more detail. In particular, we focus on the size and 
charge order of the largest cluster. Figure[2](b) shows the 
charge order parameter as a function of the size of the 
largest solid cluster along the transition path for both 
types of simulations. No significant difference is found 
between the two cases. Hence, the way the system is 
sampled does not affect the nucleation path in state A, 
confirming the hypothesis that clusters are at any time 
in quasi-equilibrium with the surrounding metastable liq- 
uid. Nothing unexpected has been observed for the study 
of state point A. 

Now we focus our attention on state point B (also at 
a supersaturation high enough as to allow nucleation to 
occur spontaneously in the course of a regular NpT sim- 
ulation). In contrast to state point A, there is a com- 



5 




5e+05 le+06 

MC cycle 



1.5e+06 




210 



FIG. 2: (a) Internal energy U versus the number of MC cy- 
cles, and (b) charge order parameter £ versus the number of 
particles in the largest solid cluster, n, for two liquid-to-solid 
trajectories, where the line connects two consecutive points 
along the trajectory. All results are obtained from two MC 
NPT simulations at state point A, where the solid line (circles) 
correspond to sampling without charge swap moves, while the 
dashed lines (squares) denote the ones with swap moves. 



petition between CsCl and disordered fee because of the 
proximity of state B to the solid-solid coexistence line. 
Figure [3] shows a different scenario as compared to the 
case of state A. Depending on the MC simulation type, 
the metastable fluid ends up forming either an ordered 
CsCl phase (with charge swap moves) or a disordered- 
fec lattice (without charge swap moves). Note that the 
energy of the disordered-fee solid (formed in the absence 
of swap moves) is higher than that of the liquid, indicat- 
ing the entropic character of the phase transition. The 
number of particles in the largest solid cluster, n, as a 
function of the number of MC cycles is shown in Fig. [3] 
(b). We have analysed the solid clusters when n starts 
to grow continuously in size (indicated by an arrow in 
FigE] (b)). Figure [3] (c), representing the charge order 
parameter as a function of the cluster size along the nu- 
cleation pathway, shows that not only a different phase 
is obtained by changing the sampling, but also that the 



nucleation path is different: if no charge swap moves are 
included in the MC simulations solid clusters have a no- 
ticeably lower charge order parameter. This is an indica- 
tion that a kinetic effect is playing an important role on 
the crystallization path. If the transition path was solely 
determined by the underlying free energy landscape, pre- 
critical clusters would have had the same charge order 
parameter regardless of the sampling scheme. 

Some of the authors have recently reported a study of 
the liquid-to-solid transition at state point C (32| , which 
is in line with the observations described here for state B. 
Since the supersaturation in state point C is lower than 
in B, a rare event simulation technique was required to 
observe the liquid-to-solid transition in computer simu- 
lations. Using the FFS technique we observed that the 
structure of the clusters along the transition path de- 
pends on the sampling scheme. If charge swap moves are 
included, the radial distribution function of the particles 
in the cluster reveals a CsCl-like structure, whereas an 
fee lattice is observed otherwise [32[ . Figure H] shows the 
charge order parameter versus the number of particles in 
the largest solid cluster for two typical transition paths 
obtained in the FFS calculations. Confirming the study 
at state point B, the path obtained with charge-swap 
moves contains clusters with higher charge order. For 
small clusters the value of £ is similar. Therefore, the 
difference between the two paths is not determined by 
the initial fluctuations towards the solid, but rather by 
the way particles are incorporated into the growing clus- 
ter. In order to grow an ordered cluster, particles have 
to be incorporated with a given charge sign at a given 
location at the cluster surface before the cluster itself re- 
dissolves. This condition is only fulfilled if swap-moves 
are included. 

In a recent paper, we have calculated the free en- 
ergy barrier associated with each type of cluster (ordered 
CsCl/disordered-fcc) via US simulations 32| . The barri- 
ers are shown in Fig. [5l We can grow either ordered or 
disordered clusters by respectively including or leaving 
out charge swaps moves in the US MC simulation. This 
allows us to calculate the free energy associated with each 
path. As Fig. [5] shows, a path containing disordered-fee 
clusters has a higher free energy. Since the US scheme is 
based on equilibrium simulations, the minimum free en- 
ergy path should be eventually reached even if the system 
is inefficiently sampled (not including charge swapping). 
In fact, for a window centred around 30 particles (the size 
beyond which the free energy associated with each type 
of path diverges), around 3 TO 5 MC cycles are required to 
complete a transition from disordered to ordered clusters. 
This "time" is much larger than the time it takes for the 
clusters to grow and follow either type of path (see Fig[3] 
(b) and (c)). Therefore, without charge-swap moves, the 
clusters keeps growing with charge disorder with no time 
to equilibrate to the lowest free energy path. 

In summary, the simulations show that both for states 
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FIG. 3: (a) Internal energy U versus number of MC cycles, 

(b) Number of particles in the largest cluster, n, versus num- 
ber of MC cycles, where the arrow indicates the point beyond 
which the liquid-to-solid trajectory has been analysed, and 

(c) Charge order parameter £ versus number of particles in 
the largest solid cluster, n, for two liquid-to-solid trajectories, 
where the line connects consecutive points along the trajec- 
tory. All results are obtained from two simulations at state 
point B, where the solid line (circles) correspond to sampling 
without charge swap moves, while the dashed lines (squares) 
denote the ones with swap moves. 
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FIG. 4: Charge order parameter £ versus number of particles 
n in the largest solid cluster for two liquid-to-solid trajecto- 
ries at state point C (both calculated by means of the FFS 
method [24l 132 ]) with (squares) and without (circles) charge 
swap moves. The line connects consecutive points along the 
trajectory. 



B and C the transition path depends on the mobility 
of the particles. A possible interpretation is that, with- 
out the charge swap moves, the fluid is not ergodic on 
the time scale of cluster growth, not allowing for a min- 
imization of the free energy of the path. The Stranski- 
Totomanow 33[ conjecture, which states that the criti- 
cal cluster structure is dictated by the lowest free energy 
barrier, is not compatible with our observations. For the 
system we present here, the mobility of the particles, and 
therefore the ability to equilibrate, not only determines 
which phase is eventually formed, but also the transition 
path. 

Since the no-charge-swap sampling is obviously more 
realistic in terms of dynamics, we expect that if an ex- 
periment were to be carried out under the same thermo- 
dynamics conditions, the system would follow the high 
free energy route of disordered-fcc clusters. We have 
confirmed by means of BD simulations -which provide 
a more realistic description of the dynamics than MC 
simulations- that a liquid quenched to state point B 
transforms into a disordered-fcc solid. 



Crystallization vs. vitrification at low temperature 
and constant volume 



So far we have presented results on the fluid-to-solid 
transition at high temperature and constant pressure. 
Let us now analyse how the fluid evolves when we quench 
the system to high u* (low temperatures) at constant 
volume (state points D-O). At the temperatures corre- 
sponding to state points A-C the solid coexists with a 
high density fluid. By contrast, it coexists with a low den- 
sity gas at the temperatures corresponding to states D-O. 
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FIG. 5: Gibbs free energy AG as a function of cluster size n 
for CsCl ordered clusters (lower barrier) and for disordered- 
fcc clusters (higher barrier) [33] . A typical configuration of 
each type of cluster (disordered-fcc or ordered-CsCl) is shown 
in the figure. Green and red particles are oppositely charged. 



FIG. 6: Fraction of solid-like particles as a function of time 
along the course of a BD NVT for states L-0 (see Fig[TJ . The 
starting configuration is a homogeneous fluid. 

100 



According to the equilibrium phase diagram, a BD NVT 
simulation starting from a homogeneous fluid quenched 
to states D-0 should show how the system phase sepa- 
rates into a low density gas and a CsCl crystal. How- 
ever, at the end of our simulations (~ 300i*) only state 
L became crystalline. In Fig. [5] the fraction of solid-like 
particles (see "Methods and technical details" section for 
the criterion to define a particle as "solid-like") is plot- 
ted against time for the quenches at packing fraction 0.5 
(states L-O). It is interesting to note that crystallization 
is only observed at state point (L). This is somewhat 
counter-intuitive, since one expects faster crystallization 
for deeper quenches (higher thermodynamic driving force 
for crystallization). Hence, it must be the slowing down 
of particles' dynamics that prevents the system from crys- 
tallizing at states M-O. We do not claim that crystalliza- 
tion will never take place at states M-0 -it will even- 
tually (as a diamond will eventually become graphite). 
We simply point out that the absence of crystallization 
within our simulation time for states M-0 is not because 
of a high nucleation barrier, but because particles can 
not move enough to rearrange into a crystalline lattice. 
Exp eriments [3], theory 16, H UHl, and simulations 
34 " 



35 



36j have confirmed the existence of an attractive 
glass for short-ranged attractive particles in a similar re- 
gion of the phase diagram. 

At lower packing fractions (cf> = 0.1 and 0.3, states 
D-K) there is no evidence of crystallization even after 
~ lOOOi*. Fig. [7] shows that the mean square displace- 
ment calculated for the structures obtained 950 t* after 
the quench develops a clear inflection point when increas- 
ing it*, which is a footprint of sub-diffusive dynamics. 
The referred slowing down of particles' dynamics could be 
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FIG. 7: Mean square displacement ((r(t) — r(0)) 2 } for the 
structures formed 950 t* after the quench to state points D-G 
(see Fig[l|. The dashed line -just for visual reference- has 
slope 1. 



responsible for the absence of crystallization for quenches 
as deep as those to states G, K or O. 

Fig. [8] shows snapshots of the system at state points 
G and K 500i* after the quench. We remind the reader 
that the initial configuration of our simulations is a ho- 
mogeneous fluid phase suddenly quenched to states G 
and K. Both configurations can be described as percolat- 
ing networks of high density amorphous branches, which 
corresponds to the typical description of a gel. For a 
percolating network of dense amorphous branches to be 
a gel, it has to possess a non zero yield stress. In a 
recent paper some of the authors have shown experimen- 
tally that oppositely charged colloids can indeed form 



gels 21|. In order to do so we imaged, by means of a 
confocal microscopy, a colloidal network while imposing 
a solvent flow through the sample. The network resisting 
the drag force exerted by the solvent demonstrates the 
ability of oppositely charged particles to form gels. In 
the same study, we demonstrated that the morphology of 
the gel-like configurations (Fig. [8]) is due to an interrup- 
tion of the metastable gas-liquid coexistence region [21 1 . 
In the equilibrium phase diagram, Fig. [IJ there is a gas- 
liquid metastable coexistence buried underneath the gas- 
solid coexistence region. When a homogeneous fluid is 
quenched to states D-K a gas-liquid spinodal decompo- 
sition starts immediately. As soon as particles gather 
in high density regions, the spinodal coarsening process 
slows down because the local dynamics of the particles 
gets slow (due, in turn, to the high energy at contact and 
the short interaction range). The restricted mobility of 
the particles not only slows down the spinodal coarsen- 
ing (giving rise to gel- like networks), but also prevents 
crystallization. 

Finally, we analyse the influence of the range of the 
potential on the ability to crystallize. This information 
can be helpful for experimentalists, for making good col- 
loidal crystals depends sensitively on numerous experi- 
mental variables [H, [37J . In order to assess the role of the 
interaction range on crystallization we have performed 
BD NVT simulations at the same packing fraction and 
u* as states L-0 but imposing a longer screening length 
(kct = 2 instead of 6). Starting from a homogeneous 
configuration we monitor the fraction of crystalline par- 
ticles over time (see Fig. |9]). Whereas for na — 6 we 
only observed crystallization at u* — 6.5 (Fig. [6]), now 
the system readily becomes solid for u* — 6.5, u* = 8 
and u* — 10. Therefore, for the same contact energy, 
increasing the interaction range favours crystallization. 
This can be intuitively understood: particles in a dense 
amorphous phase need to rearrange in order to form crys- 
talline domains; given that the attraction force is higher 
the shorter the interaction range (for a given u*), re- 
arrangements are easier for longer ranged interactions. 
From an experimental point of view increasing the in- 
teraction range in a system of charged particles implies 
decreasing the amount of salt in the medium. Fig. [9] also 
shows a manifestation of glassy particle dynamics (in the 
same way as Fig. [6] does for k<t — 6). Again, on a ther- 
modynamic basis, one would expect that increasing the 
contact energy favours the formation of the solid phase. 
On the contrary, the system becomes glassy for high con- 
tact energies and crystallization gets hindered. For the 
quench at u* — 15, a crystal cluster of ~ 250 particles 
forms at the initial stage of the simulation. The fact that 
this crystalline domain can not propagate throughout the 
simulation box is an indication that the system has fallen 
out of equilibrium. 

As shown above, at low <f> the system forms a percolat- 
ing network of dense amorphous branches when quenched 





FIG. 8: Snapshot of the structures obtained 500 t* after the 
quench to states G (a) and K (b). Green and red particles 
are oppositely charged. 
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FIG. 9: Same as Fig[6] but with a longer interaction range: 
no = 2 instead of 6. 



to a metastable gas-liquid coexistence point. In view 
of our observations for kct = 2 it is reasonable to ex- 
pect crystallization within such amorphous branches at 
a larger screening lengths. In that case, a percolating 
network of crystalline branches would result. Such a 
structure has already been observed for attractive par- 
ticles [H, [111 H(|. Some of us have recently studied 
the case of kct = 2 at packing fraction 0.1 for contact 
energies ranging from u* = 12.5 to u* = 30 [2l[. At 
u* = 20 we found that about 15% of the system was 
crystalline; for the rest of the investigated states the sys- 
tem remained amorphous by the end of our simulations. 
Nevertheless, we show in our experiments (see experi- 
mental details in the "Methods and technical details" sec- 
tion) that it is possible to obtain a percolating network 
of crystalline branches for oppositely charged particles 
(Figllip. In the experiments, we observed a percolating 
network of amorphous branches (Fig. I10|) minutes after 
putting in contact both low and high-salt concentration 
suspensions. These gels were observed in the high-salt- 
concentration half of the sample (particles are not op- 
positely charged in a free-of-salt environment 0). Two 
days later, when the salt gradient had smoothened, we 
observed that the branches had coarsened and (partially) 
cristallized (Fig. QT] (a) and [11] (b)). Hence, the forma- 
tion of percolating crystalline structures is due to local 
crystallization in the course of a metastable spinodal de- 
composition (by analogy with an amorphous-branches' 
gel, whose formation can be described as vitrification in 
the course of a metastable spinodal decomposition) . Fur- 
ther experimental work needs to be done to more accu- 
rately determine how far the time scale of gelation and 
crystallization were actually apart, as this was not inves- 
tigated yet. In contrast to the experiments shown here, in 
a recent work, we did not observe crystallization in the 
branches of oppositely-charged-particles' gels [21,]. An 



FIG. 10: A confocal microscopy image of a gel-like struc- 
ture with amorphous branches readily formed after vigor- 
ously shaking the sample with a suspension of the two particle 
species. Green and red particles are oppositely charged. 



important difference between the two experimental sys- 
tems was the amount of salt present in the medium. In 
the current experiments there is initially a salt gradient 
along the sample (0.47/zM — 0/xM). Hence, the average 
salt concentration was ~ 0.23/iM, whereas it was lfiM 
for the experiments reported in [21(. The particles' size 
was 2.5fim for the present experiments (pr) and 2fim 
for the former ones (for). Using the Derjaguin-Landau- 
Verwey-Overbeek (DLVO) [HI, [42| approximation for k 
(k = y/STrXsPsait) we can estimate the ratio between 
both kct: KO~f or /KO~p r ~ 1.7. Therefore, the interaction 
range was larger for the present experiments. This is in 
agreement with our simulations, where we observe that 
crystallization is easier for kct = 2 than for kct = 6. Sum- 
marizing, both simulations and experiments suggest that 
colloidal crystals form more readily for systems with long 
screening lengths. However, more careful experiments 
need to be carried out, since in our experimental system 
the charges of the colloidal particles are coupled to the 
concentration of salt in the medium. 



CONCLUSIONS 

We present a study of the evolution of a fluid of oppo- 
sitely charged particles after being quenched to a state 
point of the equilibrium phase diagram where either a 
solid or a solid coexisting with a gas is the stable phase. 

First of all, we have studied, by means of Monte Carlo 
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FIG. 11: (a) A confocal microscopy image of a gel-like struc- 
ture with crystalline branches observed two days after mixing, 
(b) A zoomed-in confocal image. Crystalline regions are ob- 
served in contact with amorphous ones. 

simulations at constant temperature and pressure, crys- 
tal nucleation for quenches to state points where the 
solid is the most stable phase. Interestingly, the crys- 
tallization path taken by a metastable liquid close to an 
ordered-solid/disordered-solid coexistence line depends 
on the simulations' sampling scheme. If a "fast" sampling 
scheme is used (including charge-swap moves), the nucle- 
ation path is a succession of ordered-solid clusters, while 
if a "slow" sampling is performed (not including charge- 



swap moves), the nucleation path contains disordered- 
solid clusters. By means of Umbrella Sampling calcula- 
tions we have measured the free energy associated with 
each type of path: The free energy of the ordered-clusters 
path is lower than that of the disordered clusters. We in- 
terpret our results to indicate a lack of ergodicity of the 
slowly sampled fluid on the time scale of crystal growth. 
Our results contradict the Stranski-Totomanow conjec- 
ture, which states that the transition path is determined 
by the minimum free-energy barrier for nucleation. For 
the case presented here, the mobility of the particles plays 
a major role in the selection of the nucleation path. We 
argue that if an experiment were to be carried out under 
the same thermodynamic conditions, the system would 
follow the higher free energy route of disordered clusters. 
If the same study is repeated at a state point far away 
from the ordered-solid/disordered-solid coexistence line, 
the sampling scheme does not affect the transition path. 



Secondly, we have studied the evolution of a fluid 
quenched to state points of the equilibrium phase dia- 
gram where a gas and a solid coexist. In order to do 
so, we have carried out Brownian Dynamics simulations 
at constant volume and temperature. We observe clear 
signs of slow dynamics. For instance, the system does 
not crystallize faster for deeper quenches -which is what 
would be expected if crystallization is controlled by ther- 
modynamics (nucleation free-energy barrier). Moreover, 
the mean square displacement reveals sub-diffusive dy- 
namics: a clear inflection point develops with the quench 
depth. The morphology of the system is very much influ- 
enced by the presence of a metastable gas-liquid spinodal: 
Percolating networks of dense glassy branches (gels) are 
obtained as a consequence of the arrest of an ongoing 
spinodal decomposition. Finally, we have studied the in- 
fluence of the interaction range on crystallization. In- 
creasing the interaction range, while keeping the contact 
energy and the density constant, favors crystallization 
over vitrification. We believe this is the case because 
particles can rearrange more easily the longer the interac- 
tion range is. Nevertheless, the slow-dynamics footprint 
of crystallization being delayed by an increased quench 
depth, is also present when the interaction range is in- 
creased. Crystallization in the experiments might be due, 
in agreement with our computer simulations, to an in- 
crease of screening length as the salt concentration de- 
creases. However, this has to be further investigated, 
since the amount of charge on the particles is also corre- 
lated with the salt concentration in the system. Our work 
illustrates that either amorphous or crystalline-branched 
networks can be seen as an interrupted ongoing spinodal 
decomposition. In the former case, the interruption is 
due to vitrification, whereas in the latter it is due to 
crystallization. 
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